import geopandas as gpd
import pandas as pd
from pyproj import Proj, transform
from shapely.geometry import Polygon, Point, box, MultiPolygon
import requests
import matplotlib.pyplot as plt
import numpy as np
import warnings
import csv
warnings.filterwarnings("ignore")

## 已解决多polygon问题：
class Gird_sam(object):
    def __init__(self, filepath, filename, show=False, spacing=25*1000):
        self.filepath = filepath
        self.filename = filename
        self.spacing = spacing # in meters,网格点大小
        self.show = show

    def get_points(self):
        if self.filepath.endswith('.csv'):
            address = pd.read_csv(self.filepath, encoding='gbk')
        else:
            address = pd.read_excel(self.filepath)
        csvFile = open(self.filename, 'w', encoding='gbk', newline='')  # a表示续写同名文件
        csvWriter = csv.writer(csvFile)
        csvWriter.writerow(['省', '市', '区县', 'adcode', '经度', '纬度','面积'])  #

        for i in range(len(address)):
            row = [0, 0, 0, 0, 0, 0, 0]
            adcode = str(address['adcode'][i])
            # print('adcode:', adcode)
            # https://geo.datav.aliyun.com/areas_v3/bound/530100.json
            url = 'https://geo.datav.aliyun.com/areas_v3/bound/{}.json'.format(adcode[0:6])
            # print(requests.get(url))
            name = str(address['市'][i] + address['区县'][i])  # .encode('utf8')
            print('name', name)

            try:
                res = requests.get(url).json()

                center = res['features'][0]['properties']['center']
                coordinates = res['features'][0]['geometry']['coordinates']
                points_res = self.grided_points(center, coordinates, name)
                points = points_res[0]
                areas = points_res[1]
                print(points,areas)
                row[6] = areas
                for p in points:
                    row[0] = address['省'][i]
                    row[1] = address['市'][i]
                    row[2] = address['区县'][i]
                    row[3] = adcode
                    row[4] = p[0]
                    row[5] = p[1]
                    csvWriter.writerow(row)
                    row[6]=None # 为了只写一次面积
            except:
                print(f'不存在该{name}地区的轮廓json文件')
                row[0] = address['省'][i]
                row[1] = address['市'][i]
                row[2] = address['区县'][i]
                row[3] = adcode
                row[4] = None
                row[5] = None
                row[6] = None
                csvWriter.writerow(row)

    def grided_points(self, center, coordinates, name):

        polygons = []
        area_all = 0
        gdf_utm_all = []
        utm_zone = None
        # 部分json文件的格式改变：进行修正：
        if len(coordinates)==1 and np.array(coordinates).ndim == 3:
            coordinates = [coordinates]
        for coords in coordinates:
            # todo: 没验证：
            if np.array(coords).ndim == 2:
                coords = [coords]
            polygon = Polygon(coords[0])

            # 创建 GeoDataFrame -- DataFrame that contains geometry in columns
            # crs是坐标系, EPSG:4326 就是 WGS84 的代码，WGS84是最流行的地理坐标代码（地理坐标是经纬度+海拔）
            gdf = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[polygon])

            if not utm_zone:
            # 转换多边形到投影坐标系，这里使用 UTM
                utm_zone = int((coords[0][0][0] + 180) // 6) + 1  # 计算UTM带号, 使用经度计算时区, 有的区可能跨多个，忽略误差，使用一个
            '''
            参数解释：
            * +proj:投影名
            * +zone:UTM区域
            * +north:北半球的utm区域
            * ellps:椭球体名
            * datum:基准面名
            * units:meters(米), US survey feet(美国测量英尺),等.
            '''
            proj_utm = Proj(f"+proj=utm +zone={utm_zone} +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
            gdf_utm = gdf.to_crs(proj_utm.srs)
            # 求取面积
            areas = gdf_utm.geometry[0].area
            # 多个polygon放在一起：
            area_all += areas
            polygons.append(polygon)
            gdf_utm_all.append(gdf_utm)

        multipolygon = MultiPolygon(polygons)
        all_gdf = pd.concat(gdf_utm_all)

        # 把center点转换到投影坐标系，
        lon_c, lat_c = center[0], center[1]
        wgs84 = Proj(proj='latlong', datum='WGS84')
        lon_c_utm, lat_c_utm = transform(wgs84, proj_utm, lon_c, lat_c) # 一个区的不可能跨时区，这样粗略就可以了

        # 获取转换后的多边形的边界, 转换后的边界存在了geometry[0]里:
        min_lon, min_lat, max_lon, max_lat = gdf_utm_all[0].geometry[0].bounds
        for gdf_utm in gdf_utm_all:
            min_lon_new, min_lat_new, max_lon_new, max_lat_new = gdf_utm.geometry[0].bounds
            '''
            这里要取的是公共的bounding box：therefore: 
            经度的最小值遇到更小的即更新，经度的最大值遇到更大值即更新；
            纬度的最小值遇到更小值即更新，纬度的最大值遇到更大值即更新。
            '''
            if min_lon_new < min_lon:
                min_lon = min_lon_new
            if min_lat_new < min_lat:
                min_lat = min_lat_new
            if max_lon_new > max_lon:
                max_lon = max_lon_new
            if max_lat_new > max_lat:
                max_lat = max_lat_new

        # 创建网格点
        # 假设网格间隔为 25 千米
        # 逻辑上，从center点出发，每个地区至少应该有一个点（无论spacing是多大）: 代码上是符合逻辑的
        spacing = self.spacing
        lon_coords1 = np.arange(lon_c_utm, max_lon + spacing, spacing) # 向右走
        # print('maxx + spacing:',maxx + spacing,'x_c:',x_c)
        # print('x_coords1',x_coords1)
        lon_coords2 = np.arange(lon_c_utm, min_lon, -spacing)[1:] # 向左走
        lat_coords1 = np.arange(lat_c_utm, max_lat + spacing, spacing)
        lat_coords2 = np.arange(lat_c_utm, min_lat , -spacing)[1:]
        lon_coords = np.concatenate((lon_coords1, lon_coords2), axis=0)
        lat_coords = np.concatenate((lat_coords1, lat_coords2), axis=0)
        # print('x_coords',x_coords,x_coords1,x_coords2)
        grid_points = [Point(x, y) for x in lon_coords for y in lat_coords]
        print(grid_points)

        # 检查哪些点在多边形内
        inside_points = []
        inside_points_draw = []
        for point in grid_points:
            for gdf in gdf_utm_all:
                if gdf.geometry[0].contains(point):
                    # 转换回经纬度
                    lon, lat = transform(proj_utm, Proj('epsg:4326'), point.x, point.y)
                    inside_points_draw.append(point)
                    inside_points.append((lat, lon))
        print(inside_points_draw)

        if self.show:
            bounds = all_gdf.total_bounds
            bbox = box(*bounds)
            bbox_gdf = gpd.GeoDataFrame(geometry=[bbox], crs=all_gdf.crs)
            points_gdf = gpd.GeoDataFrame(geometry=grid_points, crs=all_gdf.crs)
            inside_gdf = gpd.GeoDataFrame(geometry=inside_points_draw, crs=all_gdf.crs)

            fig, ax = plt.subplots(figsize=(10, 10))
            # 绘制地市轮廓
            all_gdf.plot(ax=ax, color='white', edgecolor='k', aspect=1)

            # 绘制边界框
            bbox_gdf.boundary.plot(ax=ax, color='green', linewidth=2, aspect=1)

            # 绘制网格点
            points_gdf.plot(ax=ax, color='red', marker='o', markersize=50, aspect=1)
            inside_gdf.plot(ax=ax, color='blue', marker='o', markersize=50, aspect=1)
            # 绘制center
            plt.plot(lon_c_utm, lat_c_utm, marker='o', color='green')

            # 设置图形显示的边界
            ax.set_xlim([bounds[0]-100, bounds[2]+100])
            ax.set_ylim([bounds[1]-100, bounds[3]+100])

            # 添加图例
            plt.legend(['Boundary', 'Bounding Box', 'Grid Points'])
            plt.title(name,loc="left")
            plt.rcParams['font.sans-serif'] = ['SimHei']

            plt.show()

        return inside_points, area_all

if __name__=='__main__':
    filepath = r'..\file\1.csv'
    # filename = r'..\file\res\1.csv'
    # filepath = '../file/修正的县级和adcode/1.csv'
    filename = '../file/修正后的网格化结果/test.csv' # 存储的文件名
    sites = Gird_sam(filepath=filepath, filename=filename, spacing=25 * 1000, show=True).get_points()

    # import os
    # path = r'D:\work_files\code\grided_district\file\修正的县级和adcode'
    # filelist= os.listdir(path)
    # for file in filelist:
    #     city = file.split('.')[0]
    #     filepath = os.path.join(path,file) # 输入的文件路径
    #     filename = f'../file/修正后的网格化结果/{city}网格化结果.csv' # 输出的文件路径
    #     sites = Gird_sam(filepath=filepath,filename=filename,spacing=25*1000, show=False).get_points()